/*  sdfrsqrt.cpp  by K.Tsuru  */
// function ID 3006     DARDIX only
/***********************************************************************
SDouble class
It provides 1/sqrt(x) using Newton's method i.e. by solving the equation
x - 1/(y^2) = 0.
[Algorism]
Taking w=x*100^p the value of 'p' is chosen to satisfy a condition 0.01<w<1.
The initial value is taken as y(0)=1/sqrt(w)(calculated by double).
It repeats
y(n+1)=y(n)+dy where  dy=0.5*y(n)(1-wy(n)^2).
Here 'y' has a value y = 1/(sqrt(x*100^p))= 1/(sqrt(x)*10^p), then considering
1/sqrt(x) = y*10^p, it multiplies 'y' by 10^p at last.
If the doubling effective figures method is used it becomes a few times faster.
It is very efficient as
  y(n) = 0.abcd ..... efgh ijkl
  dy   = 0.0000 ..... 0000 EFGH IJKL ....
the non-zero part of 'dy' overlaps with that of y(n) by a few figures.
At last it adds a correction such as
   delta = one - w*(y*y);
   if(!delta.IsHidden()) y = y + y*DsDiv(delta, 2);

**2017/08/06**** Sqrt(0.231106....) in 100000000(digits) ****** use SinTabele ******
2nd convergence
100000000  53.844(sec)  0.0 ????
          125.483     -2.0 e-100000032
          125.527       0.0
4th convergence
100000000 136.265(sec)  0.0 slower
          251.858      1.0 e-100000032
***********************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
//Multiplying by 4^n it reduces as 0.25<w<1, the precision is nearly the same.
#define reduceGT025 0   //using 0.25 < w < 1.0
#define USES_2ND_CONV_RSQRT 1 // 2nd convergence=======================
static const char* const func = "RecSqrt";
SDouble RecSqrt(const SDouble& x){
    if(x.Type() != x.REAL) x.SetError(x.RADIX_ERR, func, 3006);
    RealSize C;
    uint max_sz = x.MaxSize();
    SDouble w(x);
    double w0;

    if( x.Sign() == 0 ) x.SetError(x.DIVIDED_BY_ZERO, func, 3006);
    else if( x.Sign() < 0 ) x.SetError(x.DOMAIN_ERR, func, 3006);
    if(x.IsOne()) return ONE;   // x = 1.0

    int pow10 = w.NetRdxExp()*DFIGURES/2;   // DFIGURES is an even number.
    w.MultPowRdx( -w.NetRdxExp());  //Here 0.0001(=1/DRADIX) <= w < 1.0.

/* setting the initial value 'y0'*/
    w0 = doubleD(w);    // 0.0001 <= w0 < 1.0
#ifndef NDEBUG
assert( (w0 > 0.0) && (w0 <= 1.0 + DBL_EPSILON) );
#endif
    ulong m = 1uL;
    if(w0 < 0.01){      // 0.01 <= w0 < 1.0
        w0 *= 100.0;    m *= 100uL;     pow10--;
    }
#if reduceGT025
    ulong mpow2 = 1uL;
    while(w0 < 0.25){
        w0 *= 4.0;      m *= 4uL;   mpow2 *= 2uL;
    }
#endif
    if(m > 1uL) w = DsMult(w, m);
//It enters into the irreducible size mode.
    SDouble y(x.REAL, max_sz), delta(x.REAL, max_sz);
    y = 1./sqrt(w0);    // 'y' has about DOUBLE_FIG(=15) figures.
    /****************************/
// c:counter, itrmax:upper limit of iteration
    int c = 0, itrmax = howpow2( (DFIGURES*max_sz)/DOUBLE_FIG+1 ) +6;
    uint ef = (DOUBLE_FIG*2u)/DFIGURES, fig = C.EffFigures();
    bool fullPrec = false;  //calclation is done in full precision or not
    delta = w0 - w;
    if( delta.Sign() ){
    // If w0 is very close to double value, higher precision is necessary.
        ef = max( ef, 2u*(uint)abs( delta.NetRdxExp()) );
    }
//It takes into the fixed point mode.
    y.FixedPoint(y.RdxExp());
    if(ef > fig) ef = fig;
    do{
        if((ef = C.SetEffFig(ef)) >= fig) fullPrec = true;
#if USES_2ND_CONV_RSQRT // 2nd convergence=======================
    // faster than dev = w*y*y in DDMult()
        delta = y*y;
        delta *= w;
        delta = ONE - delta;
#ifndef NDEBUG
if(delta.FigureSize() != y.FigureSize()){
    fprintf(stderr, "%u %u\n", delta.FigureSize(), y.FigureSize());
}
assert(delta.FigureSize() == y.FigureSize());
#endif
        delta = DDiv2(delta);// delta /= 2
        delta = y*delta;
        y += delta;             //Non-zero figures overlap by a few figures.
        c++;
        ef *= 2;
        if(ef >= fig) ef = fig;
        C.SetEffFig(0);
    } while( ( !delta.IsHidden() && (c < itrmax) ) || !fullPrec );
#else // 4th convergence=======================
//cout << "  Using 4th convergence " << endl;
    SDouble h = ONE - w*(y*y);
        delta = DDiv2(h) + 0.375*(h*h);
        y *= ONE + delta;
        c++;
        ef *= 4;
        if(ef >= fig) ef = fig;
        C.SetEffFig(0);
//cout << "c= " << c << "  delta.NetRdxExp() =" << delta.NetRdxExp() << endl;
    } while( ( !delta.IsMLT(ONE) && (c < itrmax) ) || !fullPrec );
#endif
    //The condition "!delta.IsHidden()" corresponds to see the relative error.
    //The condition "c < itrmax" is not necessary but to avoid the endless repeat.
    if(c == itrmax) y.SetError(y.FATAL, func, 3006);
    y.iterationCount = c;
    y.PointFree();
    y.Reform(3006);
/*
In order to avoid "Verify" error which is rare, it adds a correction.
It evaluates d = 1-w*y'^2 where y' is an approximate value of 1/sqrt(w).
When d != 0
 y' = sqrt((1-d)/w) = y*sqrt(1-d)  ( y = 1/sqrt(w) )
 y' ~= y*(1-d/2),  y ~= y'*(1+d/2).
Then the correction is given by y'*d/2.
*/
#if 1
    if(y.PreferSpeed() == OFF){
        delta = ONE - w*(y*y);
        delta = DsDiv(delta, 2);    // delta /= 2  "delta" has about one figure.
     // delta = (1 - w*y^2)/2
        if(delta.Sign()) y += y*delta;
    }
#endif
    if(y.Verify()){     // check y = 1/sqrt(w) i.e. 1-w*y^2 << 1.0 ?
        delta = ONE - w*(y*y);
        if(!delta.IsMLT(ONE)){
            y.SetError(y.VERIFY, func, 3006);
        }
    }

    if(pow10) y.MultPow10(-pow10);
#if reduceGT025
    if(mpow2 > 1uL) return DsMult(y, mpow2);
#endif
    return y;
}